SciPy - scientific library in `python`

Adapted from the notebooks by

Introduction

SciPy builds upon NumPy.

Among others, SciPy deals with:

Integration (scipy.integrate)
Optimization (scipy.optimize)
Interpolation (scipy.interpolate)
Fourier Transform (scipy.fftpack)
Signal Processing (scipy.signal)
Linear Algebra (scipy.linalg)
Sparse matrices (scipy.sparse)
Statistics (scipy.stats)
Image processing (scipy.ndimage)
IO (input/output) (scipy.io)

%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display

Animation display with `python`

Animation with `matplotlib`

You can use FuncAnimation to animate a sequence of images:

%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")


def init():
    ax.set_xlim(0, 2 * np.pi)
    ax.set_ylim(-1, 1)
    return (ln,)


def update(frame):
    xdata.append(frame)
    ydata.append(np.sin(frame))
    ln.set_data(xdata, ydata)
    return (ln,)

ani = FuncAnimation(
    fig,
    update,
    interval=50,
    frames=np.linspace(0, 2 * np.pi, 100),
    init_func=init,
    blit=True,
)

display(HTML(ani.to_html5_video()))

Another version exists (using JavaScript):

%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "bo")

ani = FuncAnimation(
    fig,
    update,
    interval=50,
    frames=np.linspace(0, 2 * np.pi, 100),
    init_func=init,
    blit=True,
)

display(HTML(ani.to_jshtml()))

More information: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/

Animation with `plotly`

matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:

import plotly.express as px
from plotly.offline import plot

df = px.data.gapminder()
fig = px.scatter(
    df,
    x="gdpPercap",
    y="lifeExp",
    animation_frame="year",
    animation_group="country",
    size="pop",
    color="continent",
    hover_name="country",
    log_x=True,
    size_max=55,
    range_x=[100, 100000],
    range_y=[25, 90],
)
fig.show("notebook")

See more here: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/ ## Linear algebra

scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.

Documentation: http://docs.scipy.org/doc/scipy/reference/linalg.html

Solving linear systems:

Find x such that: A x = b for specified matrix A and vector b.

 A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
# A = np.random.randn(500, 500)
# b = np.ones((500, 1))
b = np.array([[1, 2, 3]], dtype=np.float64).T
# b = np.array([[1, 2, 3]], dtype=np.float64).T

print(A)
print(b)

x = linalg.solve(A, b)
print(x)
print(x.shape)
print(b.shape)

[[ 1.  0.  3.]
 [ 4.  5. 12.]
 [ 7.  8.  9.]]
[[1.]
 [2.]
 [3.]]
[[ 0.8       ]
 [-0.4       ]
 [ 0.06666667]]
(3, 1)
(3, 1)

Check the result at a given precision (different from ==)

np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)

True

Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!

Eigen values / vectors

A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:

A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)

# V = [v_1, ... v_n] , columns = eigen vectors
# A = V diag(s_1,...,s_n) V^T

np.allclose(A, evecs @ np.diag(evals) @ evecs.T)

[ 2.0582851 +0.j -2.66193205+0.j -1.8290659 +0.j] 
 ------
 [[ 0.70513863  0.70663603  0.05869437]
 [ 0.41391987 -0.47742165  0.77507349]
 [-0.57571682  0.5222395   0.62913914]]

True

EXERCISE : Eigen values/vectors

Verify numerically that the output from linalg.eig are indeed approximately eigen values and eigen vectors of matrix A above.

Hint: use https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html

If A is symmetric you should use eigvalsh (H for Hermitian) instead: this is more robust and leverages the structures (you know they are real!)

Matrix operations

linalg.trace(A) # trace
linalg.det(A) # determinant
linalg.inv(A) # Inverse, consider NEVER using it though :)

Norms

print(linalg.norm(A, ord="fro"))  # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))

3.829869697259075
3.829869697259075
2.661932052341721
[1.8290659  2.0582851  2.66193205]
3.613122169775675

EXERCISE : Norms computation

Check numerically what is the instruction linalg.norm(A, ord=np.inf) computing. Double check with the help and a numerical test.

A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))

1.9847869788550723

Random generation, distributions, etc.

(see more on this at https://albertcthomas.github.io/good-practices-random-number-generators/ or https://numpy.org/doc/1.18/reference/random/legacy.html#numpy.random.RandomState)

seed = 12345
rng = np.random.default_rng(seed)  # can be called without a seed
rng.random()

0.22733602246716966

Visualization of various popular distributions are available here (see widgets): https://github.com/josephsalmon/Random-Widgets

Optimization

Goal: find functions minima or maxima Documentation: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

from scipy import optimize

Finding (local!) minima

def f(x):
    return 4 * x ** 3 + (x - 2) ** 2 + x ** 4


def mf(x):
    return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)


xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()

Default solver for minimization/maximization: fmin_bfgs (see https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)

x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)


plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()

Optimization terminated successfully.
         Current function value: -3.506641
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8
Optimization terminated successfully.
         Current function value: -6.201654
         Iterations: 5
         Function evaluations: 12
         Gradient evaluations: 6
Optimization terminated successfully.
         Current function value: 2.804988
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8

EXERCISE : Bassin of attraction

Draw the points on the curves with two different colors : - orange for the points leading to find the left local minima - red for the points leading to the right local minima.

Find the zeros of a function

Find x such that f(x) = 0, with fsolve.

omega_c = 3.0

def f(omega):
    return np.tan(2 * np.pi * omega) - omega_c / omega



x = np.linspace(1e-8, 3.2, 1000)
y = f(x)

# remove vertical lines when the function flips sign
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)

optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))

my_zeros = (
    np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation

from scipy.optimize import curve_fit


def f(x, a, b, c):
    """f(x) = a exp(-bx) + c."""
    return a * np.exp(-b * x) + c


x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5)  # true signal
yn = y + 0.2 * np.random.randn(len(x))  # noisy added

plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()

(a, b, c), _ = curve_fit(f, x, yn)
print(a, b, c)

2.5880674567536204 1.1582704293784882 0.4083709002776674

Displaying

plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()

Rem: for polynomial fitting, one can directly use numpy

x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()

/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:

Polyfit may be poorly conditioned

Interpolation

from scipy.interpolate import interp1d, CubicSpline


def f(x):
    return np.sin(x)


n = np.arange(0, 10)
x = np.linspace(0, 9, 100)

y_meas = f(n) + 0.1 * np.random.randn(len(n))  # add noise
y_real = f(x)

linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)

cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)


plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()

Images

from scipy import ndimage, misc

img = misc.face()
type(img), img.dtype, img.ndim, img.shape


print(2 ** 8)  # uint8-> code sur 256 niveau.

n_1, n_2, n_3 = img.shape
np.unique(img)


# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()

/tmp/ipykernel_25871/4081896913.py:3: DeprecationWarning:

scipy.misc.face has been deprecated in SciPy v1.10.0; and will be completely removed in SciPy v1.12.0. Dataset methods have moved into the scipy.datasets module. Use scipy.datasets.face instead.

RGB decomposition

fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))

plt.tight_layout()

print(img.flags)  # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags)  # can edit now


img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show()

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : False
  ALIGNED : True
  WRITEBACKIFCOPY : False

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

Convert a color image in grayscale

plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()

Blurring (fr: floutage)

img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()

Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)

%matplotlib widget
import ipywidgets as widgets

# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()

@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
    """Remove old lines from plot and plot new one"""
    # [l.remove() for l in ax.lines]
    img_flou = ndimage.gaussian_filter(img, sigma)
    ax.imshow(img_flou, cmap=plt.cm.gray)

Changing colors in an image

import pooch
import os

url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)

img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)


fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)

plt.tight_layout()

EXERCISE : Make the Brazilian italianer

Create a version of the Brazilian flag as follows:

SciPy - scientific library in `python`

Introduction

Animation display with `python`

Animation with `matplotlib`

Animation with `plotly`

Solving linear systems:

Eigen values / vectors

EXERCISE : Eigen values/vectors

Matrix operations

Norms

EXERCISE : Norms computation

Random generation, distributions, etc.

Optimization

Finding (local!) minima

EXERCISE : Bassin of attraction

Find the zeros of a function

Parameters estimation

Interpolation

Images

RGB decomposition

Convert a color image in grayscale

Blurring (fr: floutage)

Changing colors in an image

EXERCISE : Make the Brazilian italianer

More on colors/ images:

Additional lectures

Introduction

Animation display with python

Animation with matplotlib

Animation with plotly

Solving linear systems:

Eigen values / vectors

EXERCISE : Eigen values/vectors

Matrix operations

Norms

EXERCISE : Norms computation

Random generation, distributions, etc.

Optimization

Finding (local!) minima

EXERCISE : Bassin of attraction

Find the zeros of a function

Parameters estimation

Interpolation

Images

RGB decomposition

Convert a color image in grayscale

Blurring (fr: floutage)

Changing colors in an image

EXERCISE : Make the Brazilian italianer

More on colors/ images:

Additional lectures

Animation display with `python`

Animation with `matplotlib`

Animation with `plotly`